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Abstract 

A recipe for creating materials with a desired refraction coefficient is 
implemented numerically. The following assumptions are used: 

( m = h(x m )/a K , d = 0(a (2 ~ K)/3 ), M = 0(l/a 2 ~ K ), k G (0, 1), 

where £ OT and x jfi are the boundary impedance and center of the m-th 
ball, respectively, h(x) £ C(D), lmh(x) < 0, M is the number of small 
balls embedded in the cube D, a is the radius of the small balls and d is 
the distance between the neighboring balls. 

An error estimate is given for the approximate solution of the many- 
body scattering problem in the case of small scatterers. This result is used 
for the estimate of the minimal number of small particles to be embedded 
in a given domain D in order to get a material whose refraction coeffi- 
cient approximates the desired one with the relative error not exceeding 
a desired small quantity. 

MSC: 65R20, 65Z05, HQ 10 

Key words: many-body wave scattering problem, metamaterials, refraction co- 
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1 Introduction 



A theory of wave scattering by many small bodies embedded in a bounded 
domain D filled with a material with known refraction coefficient was developed 
in PQ-[1]. It was assumed in Q] that 

d = 0(a 1 / 3 ), M = 0(ar 1 ), — ^— - = Cm u M on S m , 1 < m < M, 

ov 

where a is the characteristic size of the small particles, d is the distance between 
two neighboring particles, M is the total number of the embedded particles, S m 
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is the boundary of m-th particle D m , v is the unit normal to S m directed out 
of D m , and Cm = h m /a, where h m , Imh m < 0, 1 < m < M, are constants 
independent of a. 

Let us assume that D is filled with a material with known refraction co- 
efficient ng(a;), Imn^i) > 0, Hq(x) = 1 in D' := R 2 \ D, n^(x) is Riemann 
integrable. The governing equation is 

L u := [A + k 2 nl]u = 0, in K 3 , (I) 



u.q = cxp(ikx • a) + Vq (2) 

where k is the wave number, a G S 2 is the direction of the incident plane wave, 
S 2 is the unit sphere in R 3 , and vo is the scattered field satisfying the radiation 
condition 

lim r [ ~g— ~ ikvo ] =0 r := \x\ — > oo, (3) 



r— f oo 



and the limit is attained uniformly with respect to the directions x° :— x/r. 

Let n 2 (x) be a desired refraction coefficient in D. We assume that n 2 (x) is 
Riemann integrable, Im n 2 (x) > 0, n 2 (x) = 1 in D' . Our objective is to create 
materials with the refraction coefficient n 2 (x) in D by embedding into D many 
small non-intersecting balls B m , 1 < m < M, of radius a, centered at some 
points x m E D. If one embeds M small particles B m in the bounded domain 
D, then the scattering problem consists of finding the solution to the following 
problem: 

Lqu m := [A + k 2 n 2 ]u M (x) =0 i£l 3 \ U^ =1 B m , (4) 

Bum 
di 



= ( m "M on S m := dB mi 1 < m < M, (5) 



u M = u + v M , (6) 

where uo solves problem ([J)-©, and vm satisfies the radiation condition. 
The following theorem is proved in [3] under the assumptions 

( m = h(x m )/a K , d = 0(a (2 - K >/ 3 ), M = O {l/a 2 - K ) , ke(0,1), (7) 

where h(x) is a continuous function in D, Imh < 0, and k € (0, 1) is a parameter 
one can choose as one wishes. Below it is always assumed that conditions ([7]) 
hold. 

Theorem 1.1 ([4 j). Assume that conditions ([7]) are satisfied, and D m is a ball 
of radius a centered at a point x m . Let h(x) in ([7J be an arbitrary continuous 
function in D, Imh(x) < 0, A p c D be any subdomain of D, andJ\f(A p ) be the 
number of particles in A p , 

7V(A p ) = --L / N(x)dx[l+o{l)}, a^0, (8) 
a Ja„ 
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where N(x) > is a given continuous function in D. Then 

lim \\u e (x) -u(x)\\c(D) =0, (9) 

a— fO 

ui/iere 

M 

u e (x) := u (x) - Air'^G{x,x : j)h{x j )u e {x j )a 2 ~ K [l + o(l)], a ->• 0, (10) 

where mirij \x — Xj\ > a, and G(x,y) is the Green function of the operator Lq in 
R 3 , G(x,y) satisfies the radiation condition. The numbers u e (xj), I < j < M, 
are found from the linear algebraic system: 

M 

u e (x m ) — u (x m )-4ir ^2 G(x m ,x J )h(x J )u e (x : j)a' 2 ^ K , m = 1,2, . . . , M, 

3=1 ,j#m 

(11) 

which is uniquely solvable for all sufficiently large M. The function u{x) = 
lim a _j.o Ue (x) solves the following limiting equation: 



u(x) = u (x)- G(x,y)p(y)u(y)dy, (12) 

JD 

where uq satisfies equations ([l])-([3]), 

p(x) := 4irN(x)h(x), (13) 

n 2 {x) := 1 - k- 2 q(x), (14) 

q(x) := q (x) +p(x), q (x) := k 2 - k 2 nl(x), (15) 
and TIq(x) is the coefficient in ([1]). 

In [4] a recipe for creating material with a desired refraction coefficient is 
formulated. 

The goal of this paper is to implement numerically the recipe for creating 
materials with a desired refraction coefficient in a given domain D by embedding 
in D many small particles with prescribed physical properties. These particles 
are balls of radius a, centered at the points x m € D, and their physical properties 
are described by the boundary impedances £ m = h(x m )/a K . A formula for 
embedding the small balls in D is given in Section 2. 

We give an estimate for the error in the refraction coefficient of the medium 
obtained by embedding finitely many (M < oo) small particles, compared with 
the refraction coefficient of the limiting medium (M — > oo). This is important 
because in practice one cannot go to the limit M — > oo, i.e., a — > 0, and one 
has to know the maximal a (i.e., minimal M) such that the corresponding to 
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this a refraction coefficient differs from the desired refraction coefficient by not 
more than a given small quantity. In Section 3 we give an algorithm for finding 
the minimal number M of the embedded small balls which generate a material 
whose refraction coefficient differs from a desired one by not more than a desired 
small quantity. In Section 4 some numerical experiments are described. 



2 Embedding small balls into a cube 

In this section we give a formula for distributing small balls in a cube in such a 
way that the second and third restrictions (J7J are satisfied. 

Without loss of generality let us assume that the domain D is the unit cube: 

D := [0,1] x [0,1] x [0,1]. (16) 

Let 

D = U^ATg, n G N, A, n A 3 =0 for i £ j, (17) 

where N is the set of positive integers, X is the closure of the set X, and A q , 
q = 1, 2, . . . , n 3 , are cubes of side length 1/n. 

Definition: We say that D has property Q n if each small cube A q contains 
a ball of radius a n , < a n < 1/n, centered at the centroid of the cube A q , and 
the following condition holds 

d n := min dist(B an (x g ), B an ( Xj )) = ia\ 2 -^' 3 , (18) 
where x q is the centroid of the cube A q , q = 1, 2, . . . , n 3 , 

B a {x) := feel 3 | \y-x\< a}, (19) 

and 7 > is a constant which is not too small (see formula (|25[) ). 
From dXTj) and flTSJl one gets 

d n = l n - 2a n = 7 ai 2 ~ K)/3 , In ■= 1/n. (20) 
Since l n — 1/n, the quantity a n solves the equation 

7a (2-«)/3 + 2a - 1/n = 0. (21) 

The function f(a) := ja( 2 ~ K ^ 3 + 2a is strictly growing on [0, oo). Thus, the 
solution to equation (|21[) exists, is unique, and can be calculated numerically, 
for example, by the bisection method. 

However, it is easy to derive an analytic asymptotic formula for a n as n — > oo. 
This formula is simple and can be used for all n we are interested in, since these 
n are sufficiently large. 

Let us derive this asymptotic formula. Since 1/3 < (2 — k)/3 < 2/3, one has 
a < a^ 2 -^/ 3 if a < 1. Therefore, 

o„ = [l/{ ni )fl {2 ~ K \l + o(l)], asrwoc, (22) 
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is the desired asymptotic formula for the solution to (f^Tj) . Note that 
lim a n = and lim na n = 0, (23) 

as follows from (f22j) because 3/(2 — k) > 1. 

We note that a n /d n <C 1, if n » 1, because (f2T)]) yields 

On/dn = a n /(74 2 " K)/3 ) = 4 1+K)/3 /7 « 1- (24) 
Let us choose n sufficiently large so that 

7 » (Z„/2)( 1+K >/ 3 , k (£ (0, 1), Z„ = l/n, (25) 
and make the following assumption: 

Assumption A): D has property Q m p. Here D = Ug =1 fl q , fljDili = for j ^ i, 
where each cube f2 g has side length 1 / P , and m 3 small balls are embedded in Q. q 
so that the following two conditions hold: 

1. Each cube fl q is a union of small sub-cubes Aj i? : 

n q = U^Aj, „ A„nA M =0 for i ^ j, (26) 

where Aj ig , j = l,2,...,m 3 , q — 1,2,...,P 3 , are cubes of side length 
l/(mP), ' 

2. In each sub-cube Aj ig there is a ball of radius a m p, < a m p < l/(mP), 
centered at the centroid of the sub-cube Aj jg , and the radius a m p of the 
embedded balls satisfies the relation 

1/K) - 2a mP = 7^ P K)/3 , 7 » [l/(2mP)]« K+1)/3 , (27) 

where 7 > is a fixed constant. 
Lemma 2.1. If Assumption A) holds, then 
lim Ma 2 m ; = I/7 3 , (28) 

m— foo 

where M — (mP) 3 is i/ie iota/ number of small balls embedded in the unit cube 
D, 7 > is fixed, and 

a mP = [l/(7mF)] 3/(2 " K) [l + o(l)] as m — y 00. (29) 

Proof. Relation (f2"9"| is an immediate consequence of (|2T|) . Using this relation, 
one obtains 



lim Ma^ = lim (mP) 3 a^? 

= lim (mP) 3 [l/( 7 mP)] 3 [l + o(l)] (30) 

m— voo 

= lim (l/ 7 3 )[l + (l)] = l/ 7 3 . 

m— >oo 

Lemma l2.1l is proved. □ 
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3 A recipe for creating materials with a desired 
refraction coefficient 

In this section the recipe given in [3] is used for creating materials with a desired 
refraction coefficient by embedding into D small balls so that Assumption A) 
holds. 

Step 1. Given the refraction coefficient tiq(x) of the original material in D 
and the desired refraction coefficient n 2 (x) in D, one calculates 

p(x) = k 2 [n 2 (x) - n 2 (x)} = p t (x) + ipz(x), (31) 

where 

Pi := Re p(x) and P2(x) := Imp(x). 
Choose 

N(x) = l/ 7 3 , (32) 

where 7 is the constant 7 in Assumption A). 
Step 2. Choose 

h(x) = hi(x) + ih 2 (x), (33) 

where the functions hi(x) and h,2{x) are dehned by the formulas: 

K{x) = 7 3 p i (^)/(47r) ) * = 1,2, (34) 

and the functions Pi(x) are defined in Step 1. 

Step 3. Partition D into P small cubes f2 p with side length 1/P, and embed 
m 3 small balls in each cube il p so that Assumption A) holds. 

Then 

Af(Q p ) = * / N(x)dx = P P \/tfc?-p) = 1 /[ 7j P^p K)/3 ] 3 , (35) 



where A/"(A p ) is the number of the balls embedded in the cube Q p , k £ (0, 1), 
a m p is the radius of the embedded balls, and |f2 p | is the volume of the cube £l p . 
Since 

a m p = [l/(rn7P)] 3 ^~[l + o(l)] as m -> 00, 
it follows that 

hm ^1 = 1. (36) 

2-" 1+0 

By Assumption A) the balls are situated at the distances 7a m 3 P , 7 > [1/ (mP)] 3 
Therefore, all the assumptions, made in Theorem lf.il hold. Thus, 

max I ti e (a;) — u(x)\ as M —¥ 00, (37) 



G 



where u e (x) is defined in (ITU1) and u(x) solves (JT2J) . Let us assume for simplicity 
that tiq(x) = 1, so that 



y) = g(x, y) := exp(ifc|a; - y|)/(47r|a; - 



Then 



M 

u e (x) = n (a:)-47ry^g(a;,a;j)/i(xj)u e (a:j)a^'p , |ac— > a mP , M := (mP) 3 , 

(38) 

and the limiting function 



u(x) = lim u e (x) 

M— ¥00 

solves the integral equation 
u(x) + Tu(x) = uq(x), 
where 

Tu(x) := / g{x,y)p{y)u{y)dy, 
Jd 



g(x, y) := exp(ifc|a; - y\)/{An\x -y\), 

M := (mP) 3 , h(x) = h 1 (x) +ih 2 (x), hi(x),i = 1,2, are defined in (1541 . 



(39) 

(40) 
(41) 



p(x) = k 2 [nf,{x) - n 2 (x)] = 4ir[hi(x) + i/i 2 (x)]iV(a;) = 4tt[/h(x) + i/i 2 (a;)]/r , 



and the function uo(x) in (13"8l solves the scattering problem (p}-©. 
It follows from (l37l) that 



max litfe) — u e (xi)\ — >■ as M — > 00, 

KKM 



(42) 



(43) 



where w e (x) and u(x) are defined in (|3"51) and (|39p . respectively. Here and 
throughout this paper D := UfLxDj, M := (mP) 3 , Dj (j = 1,2,...,M) are 
cubes with the side length l/(mP), Dj n -D; = for j ^ I, and Xj denotes the 
center of the cube Dj. However, since (|37|) was not proved here, let us prove 
relation (031. We denote 



( Ul \ 

V2 



\u\\oo := sup \u(x)\ and ||u||, 
16D 



■M := max ti, , v 

l<j<M 1 31 
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where C is the set of complex numbers. 

Consider the following piecewise-constant function as an approximate solu- 
tion to equation 



M 



J'=l 



M ■ 



where Uj m (J = 1,2, ... , M) are constants and 
Xj(x) ■= 



1, x £ Dj, 
0, otherwise. 



(44) 



(45) 



Substituting U(m){x) for u(x) in ([39]) and evaluating at points xi, one gets the 
following linear algebraic system (LAS) which is used to find the unknown Uj^[\ 



U(M) + Td,MU(M) = W0,M, 

where m is a discrete version of Tm , defined below, 



U2,M 



t(M) 



e C M , u ,m 



\ Um,M J 

uq(x) solves problem (l)-(3), and 



I Uq{xi) ^ 
U (x 2 ) 

\ u (x M ) J 



(46) 



(47) 



M . 

(Td,Mv)i -=y2 g(xi,y)p(y)dyvj, I = 1,2, . . . ,M, v 



''2 



(48) 



Multiplying the Z-th equation of (|46|) by Xz(a:), 2 = 1,2, ... , M, and summing 
up over I from 1 to M, one gets 



«(M)(aO = Uo,(M)(x) - T M u {M) (x), 
where i^m)^) is dchncd in (1441) . 

M 

u ,(M)(x) :=J2xj(x)uo(xj), 



and 



(49) 



(50) 



(51) 



T M u(x) :=y2xj(x) / g{xj,y)p{y)u{y)dy, M=(mP) 3 . 
i=i Jd 

It was proved in [5] that equation (|49| is equivalent to (l46l) in the sense that 
{wj,m}^i solves ((46j) if and only if function (jSJ) solves (|49|) . 
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Lemma 3.1. For all sufficiently large M equation (|46[) has a unique solution, 
and there exists a constant c\ > such that 

|| (Id.M + Td.M)" 1 1| <ci, VM>M , (52) 
where Mq > is a sufficiently large number. 

Proof. Consider the operators T and Tm as operators in the space L°° (D) with 
the sup-norm. Let ||(T — Tm)m||oo :— swp xeD \(T — Tm)u(x)\. Then 

M „ 

||(T-Tm)w||oo < max sup V / \(g(x,y) - g(x u y))p(y)u(y)\dy 
4 i6A f— r J_d, 



1 l»-^l<^p j = l J °J 



(53) 



< 0(l/(mP)). 

This implies 

||T - T M \\ = 0(l/(mF)) = 0(l/M 1/3 ) -> as M -> oo. (54) 

The operator 7 + T is known to be boundedly invertible, so ||(7 + T) _1 || < c, 
where c > is a constant. Therefore, 

I + T M = (I + T)[I + (7 + T)-\T M - T)\. (55) 
By (fS~4"|) there exists Ma such that 

||(7 + T)- 1 (Tm-T)||<c||T m -T||<<5<1, VM > Af , (56) 
where (5 > is a constant. From (|5"6")l we obtain, VM > Mo, 

|| [7 + (7 + T)-HT M - T)]-|| < rT^TTTTHTWOl * V(l - *)■ (57) 
Therefore, it follows from (|55")l that 7 + Tm is boundedly invertible and 
(7 + Tm)- 1 = [7+(7 + T)- 1 (T M -T)]- 1 (7 + T)- 1 , (58) 
so there exists a constant Co > such that 

j|(7 + T M )- 1 || <c , VM>M . (59) 

Since f4"9")) is equivalent to (|4"6")l . it follows that the homogeneous equation v + 
Td,MV = has only trivial solution for M > M , i.e., Af(Id,M + Td.Af) = {0} 
for M > Mo, where Af(A) is the nullspace of the operator A, Id,M is the 
identity operator in C M and Td,M is defined in (|4"5)l . Therefore, by the Fredholm 
alternative equation (|46"]) is solvable for M > Mo. This together with (|59|) yield 
the existence of a constant ci > such that || (Id.M+Td.M)^ 1 1 < Ci for M > M . 
Lemma [3~T1 is proved. □ 
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Define T d : C 2 {D) -> C M as follows 

M . 

{T d w)i := (Tw)(xi) = y2 g( x hy)p(y) w (y)dy, 



1 = 1,2,. ..,M, 



(60) 



3=1 



and 



MM 



/ w(a;i) ^ 
u(i 2 ) 



\ u(x M ) / 



(61) 



where T is defined in (gDJ) and u(ac) solves ([32) ■ Then it follows from (pOj) 
and (I5T1) that the following equation holds 



u M + T d u = ti . 



M ■ 



(62) 



where Uo,m is defined in (|4"7|) . Using equations and (|5^t . we derive the 
following equality: 

(Jd,M + T d: M)(U(M) — u m) = {Id,M + ?d,M)u(M) — Ud,M + T^m) u M 

= u a .M - (Id, mum + T^ m um) (63) 

= Uq,M — Um — Td,M u M = TdU — Td,M u Mi 

where u/m) an d "o,Af are defined in (|47p. 



V2 



G C 



a/ 



\ DM / 

Using relation (|63"|) , one gets 

U(m) - u M = (Id,M + Td,M)~ 1 (TdU - T^ m um)- 



(64) 



(65) 



Lemma 3.2. Let Assumption A ) hold (see Section 2 below (|25l) ). Suppose u(x) 
solves (|39p where p(x) G C 1 . TTien 



- "m||c« = 0{l/M 2/3 ) as M -> oo, 
where Um and U(m) are defined in (1611) and (|47[) . respectively. 
Proof. By (|65|) and estimate (|52t we obtain 

\\U(M) - U M \\c M < ll(^d,M + Td,M)~ || \\TdU ~ T^AfUMlIc^ 
< Ci||T d U - Td,MUM ||<C M > 



(66) 



(67) 
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where T^m and T d are denned in (|35]l and respectively. Using the identity 
— u(x.j) = u(y) — u(xj) — T>u(xj)(y — Xj) + T>u(xj)(y — Xj) and applying 
the triangle inequality, we get the estimate 



\\T d u - T d M u M \\c M = max 

KKM 



< max 

KKM 



max 

KKM 



D 
M 



Y / 9(xi,y)p{y){u{y) - u(xj))dy 
g(xuy)p(y)(u(y) - u{xi))dy 

Y / 9(xi,y)p(y)(u(y) -u(xj))dy 
1 , „Jfli 



< 2||p|| ooljuHoc max / |g(xz, y)\dy 



KKM 



max 

KKM 



<2||p||ooll« 



Y / 9(xi,y)p{y)(u(y) -u(xj))dy 

r/i/ 



max 

KKM 



M . 

Y / g(xi,y)p{y){u(y) - u(xj))dy 



2||p|| 00 ||it|| 



v / 3/(2mP) 



rcir 



max 

KKM 



< 



3|bHoo||tx| 
2{m.Pf 



M . 

Y / 9{xi,y)p{y){u{y) ~ u{x 3 ))dy 



h + h 



where 



1 1 := max 

KKM 



and 

I2 := max 

KKM 



M 



3 = 1,3^1 



M 



Y / 9(xi,y)p{y){u{y) - u(xj) ~Vu( Xj )(y ~ Xj j)dy 



(68) 



(69) 



(70) 



Let us derive an estimate for I\. Using the Taylor expansion, we get 

M „ 

h < IblU max V / \g(xi,y)\ sup |X> 2 u(sy+(l-s)a; :) )||2/-a; 3 | 2 dy. (71) 
i<z<Af .rr;,,JDi o<s<i 
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Since p <E C^D), u € C 2 {D), J D \g(x,y)\dy < oo and \y-Xj\ < for y e D 
it follows from (|7T|) that 



/i = 0(l/(mP) 2 ) = 0(1/M 2/3 ), as M 



oo. 



(72) 



Estimate of I 2 is obtained as follows. Since Xj is the center of the cube Dj, it 
follows that 

g(xi,Xj)p(xj)T>u(Xj)(y-Xj)dy = 0, j = l,2,...,M. (73) 
Therefore, using ([75)l . J2 can be rewritten as follows: 



I2 = max 

KKM 



I (9(xi,y)p{y) - g{x h Xj)p{xj))Vu{xj){y - x 3 )dy 



-D 1 



Let 

m(y) ■= g(xi,y), {gip){y) = gi(y)p(y), 1 = 1,2,. ..,M. 

Then the formulas 



\(gip)(y) - {gw){xj)\ = 



1 a 

jr,{gw){ty + {i-t)x 3 )dt 

01 



< sup \V y (gip)(ty + (1 - i)xj)| |y - Xj\ 

0<t<l 

and V y (gip)(y) = p{y)V y gi(y) + gi(y)Vp(y) , yield the following estimate: 
\/3||X>u| 



/ 2 < 



max 



< 



2mP \<1<M 

\/3||Pm||oo 

2mP 
c(k)\\Vu\\ c 



sup \V y gi(ty + (1 - t)xj)\\y - Xj\dy 

D 0<t<l 



(mP) 
c||x>m||oo 
(mP) 



max HUplloo / sup |g i (i2/+(l-i)a;j)||y-x J -|(iy 
1<(<M Jdo<ki 



max p oo 
2 Kl<» 



1 



1 



(74) 



(75) 



(76) 



D \4ir\xi - y\ 4n\xi - y\ 2 

, max ||Dp|U / — - ,^0(l/(mP) 2 ) = 0(l/M 2 / 3 ), 

- l<i<M J D 47r|x; - y| 



(77) 



where c > is a constant and c(k) is a constant depending on the wave number 
k. Here the estimates \y — Xj\ < \/3/(2mP) for ye Dj, J D 4 „i x ^_ y i0 dy < 00 for 
jH < 3, and |xj - y| < 2|x/ - s| for y e Dj, j ^ I, s = txj + (1 - t)y, t e [0, 1], 
were used. The relation (JgBJl follows from (IFfl). (|68 | . (|72|) and (|77|) . 
Lemma [3~2l is proved. □ 
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Lemma 3.3. Let the Assumption A) hold. Consider the linear algebraic system 
for the unknowns u e (xi): 



,(xi) = u (xi) - 47r 22 9{xi,Xj)h{xj)a 2 m pU e ^ 

3=1,3*1 



3 )■> 



Z = 1,2,...M, 



(78) 



where p(x) = <±iih(x)N(x) £ C 2 (D), N(x) = I/7 3 , M = (mP) 3 , and g(x,y) is 
defined in (|40[) . Then 

\\u(M) ~ u eM \\ c m = O {^L + |1 - 7 3 ^"Cp l) as M ->■ °°> (79) 
where U(m) * s defined in (I4T[) . 



/ Ue(xi) \ 

Ue(x 2 ) 



U e ,M 



(80) 



\ u e (x M ) j 
and u e (xj), j — 1,2, ... , M, sotoe system ([78 
Proof. Let us rewrite (f78| as 



,1/ 



3 (x;) = ito(^) - g(xi,Xj)p(o 



r U e (^)|^ 



(81) 



u (a;;) - (T e u e . M )i, I = 1,2,...,M, 



where p(x) = Airh(x)N{x), \Dj\ = l/(mP) 3 is the volume of the cube Dj, 
N(x) = I/7 3 and 



M 3 

(T e v)i := ^2 9{xhX j )p{x j )[-fmPaf nP K)/ ^ \Dj\vj, v = 



3=1,3*1 



( ^ \ 



-1 A/ 



(82) 



/ = 1,2, ...,M. 

Let us derive an estimate for ||u(m) — u c ^m\\c m ■ Using equations (|46[) and (|81[) . 
we obtain 

(-^2,M + Td.M){u(M) — U ej MJ — {Id,M + Td^Al)u(M) — (Id,M + T d ^ M )u e ,M 

= tio,M _ Ue , M — Td.MU e ,M = T e Ue,M ~ Td,MU e ,M , 

(83) 
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where U(m) and uo : m are defined in (|4"T1) . u e ,M, Td,M and T e are denned in ([8"0jl. 
(|48|) and (|82l) . respectively. Relation ([83]) implies 



(84) 



|M(M) — "e,M||c M < IIUd.M + T^m) 1 || ||T' e M e! M — Td,MU e ,M ||c M 



< Ci||T e M e M - T d M u 



e,M\\C M i 



where estimate ([521) was used. 

Let us derive an estimate for \\T e u ei M — Td t MU e ,M\\c M ■ Using definitions 
(|4"5|) and (|82p . and applying the triangle inequality, one gets 

|| (T e — Td y M)u e ,M \\c M 



= max 

KKM 



max 

KKM 



< max 

KKM 



max 

KKM 



max 

KKM 



E 9i{x 3 )Pa mP {x 3 )u e {x 3 )\D 3 \-Y I gi(y)p(y)dyu e (x 

/ (9l{xj)Pa mP (xj) - gi{y)p{y))dyu e {xj) -I gi(y)p(y)dyu e (xi) 

1 --jJDi Jft 



gi(y)p(y)u e (xi)dy 







+ max 

1<2<M 





M 



X! Pompfe) / imiv) - m{x 3 ))u e {xj)dy 



D, 



< max (J o (0 + J X (Z) + J 2 (Z) + J 3 (l)), 

KKM 



where cft(y) := g(x t ,y), Pa mP (x) := p(x) (imPa < £ p ^ /3 ^ 
M l ) '■= f \g(xi,y)p(y)u e (xi)\dy, 

g{xi,y){p(y) - p(xj))u e (xj)dy 





Jd, 






M 




Jx(Z) 


■= E 


/ 








M 




■MO 










and 


M 




Mi) 





{g{xi,y) - g(xi,Xj))u e (xj)dy 



(85) 

(86) 
(87) 
(88) 

(89) 
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Using the estimate \xi — y\ < v3/(2mP) for y 6 Di, one gets the following 
estimate of Jq(1): 



Mi) ^ IblloolKellc-w / 



< 



Ibll oo I be || C M 



v / 3/(2mP) 



rrfr |b||oo|be||c M 



gbUH^f =0(1/M 2/ 3) 



where B a (x) is defined in (1191) . 

Let us estimate J%(1). Using the identity 

P(y) - P( x j) = P(y) - P{xj) - T>p(xj) ■(y-x j )+ Vp(xj) ■ (y - Xj) 
in (|87p and applying the triangle inequality, one obtains 

Ml) < Jl,l + Jl,2, 

where 



M 

Jia--= E 

and 

M 

E 



D, 



g{xi,y)[p(y) -p( x 3 ) - Vp(xj) ■ {y - xj)]u e (xj)dy 



g(xi,y)T>p(xj) ■ (y - Xj)u e (xj)dy 



D, 



(90) 



(91) 
(92) 

(93) 
(94) 



To get an estimate for Ji t i, we apply the Taylor expansion of p(x) and get 



Ji,i < 



be||c« 



M . 

E / \9( x i:V)\ SU P \V 2 p{ty + (1 - t)xj)\\y- Xj\ 2 dy 

T~,,JDi 0<t<l 



< 



3||^ 2 p||oo|be||c" 



M 



E 

3=1,3^1 



\g{xi,y)\dy 



(95) 



8(mP) 2 

= 0{l/{mPf) = 0(1/M 2/3 ) asM^oo, 

where B a {x) is defined in ([T^|) . and the estimate \y — Xj\ < \/3/(2mP), y G Dj, 
was used. 

Using the identity 



g(xi,Xj)Vp(xj)(y - Xj)u e (xj)dy = 0, j = 1, 2 



....,M, 



(96) 
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one derives the following estimate of 

M 

(g(x h y) - g(x h x J ))Vp(x j ) ■ (y - Xj)u e (xj)dy 



M 



< IK||c« Y] / \(9(xi,y)-g(xi,x j ))Vp(x j )-(y-x j )\dy 

M 



< 



< 



\/3||p||oo||Me||c M 



2mP 



E / \g(xi,y) - g(xi,Xj)\dy 



oo || || C M 



M 



2mP 



V / sup \T>g(x h ty + (l-t)x j )\\y-x j \dy 



(97) 



< c(fc)lbl|oo|ke|lc M 



< 



(mP) 2 

c(fc)||p|| 00 ||M e || C Af 

(toP) 2 

c ( fc )lbl|oo||w e ||c M 



M 

E 



sup 



(mP) 2 



Dj 0<t<l 47r|x ; - iy- (1 - t)x-, | 2 

1 

sup j — i 

Oi o<t<i ^\xi - ty - (1 - i)£E,-| 

1 1 



dy 



s^cxi) v 47r l^ - y\ ^Vi - ?yp 



dy 



0(l/(mP) 2 ) = 0(1/M 2/3 ) as M — > oo, 



where c(fc) is a constant depending on the wave number k and B a (x) is defined in 
(fT9|) . Here the estimates \y — Xj\ < \/3/(2mP) for y e Dj, J D 47r \ x f_ y \$ dy < oo 
for /? < 3, and |xj — y\ < 2|xj — s| for y € D^, j ^ I, s = txj + (1 — i)y, i € [0, 1], 
were used. Applying estimates (f9"51 and (liJTl) to (|9^|) . we get 



as M — » oo. 



(98) 



Let us derive an estimate for J2(Z). From ((88)) and the definition Pa mP {x) 
p(x)('jmPa < - 2 ~ K ^ 3 ) 3 we get 



■MO < ||u e ||cM 



E [ |5(x i ,y)|b(x,)||l-( 7 mPa( 2 - K )/ 3 ) 3 |d 2/ 

<||u e || c «|b|| 0o |l-( 7 mPa( 2 -^/ 3 ) 3 | T / \g(x h y)\dy 

dy 



< lluellcMllplUll-CTmPo^-^si 



rdr) || Me || C M||p||oo|l-(7^Pa (2 - K)/3 ) 3 | 



= -|| Ue || c „|b|| 00 |l- 7 3 Ma^ ) |=0(|l- 7 3 Ma^ ) |), 



(99) 
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where B a (x) is defined in p^|) . Using the relation liniM->-oo 7 3 Ma' 2_K ' = 1, one 
gets J2(0 — > as M — > oo. Estimate of J^H) is derived as follows. Using the 
identity 



(100) 



Vg(x h Xj)(y - Xj)u e (xj)dy = 0, j = 1,2, ...M, 
one gets the following estimate: 



M 



js(o= E 



(xj) / \g(xi,y)-g(xi,Xj)-Vg(xi,Xj)-(y-Xj)]dy 



< 



|b|U 7 3 Ma^ A 



D, 



sup \T> 2 g(x h ty + (1 - *)xj)||y - ^ | 2 d?/ 



< 



< 



< 



< 



< 



3c M 
.(mP) 2 

c{k)c M 
(mP) 2 

c(k)c M 
(mP) 2 

c(k)c M 
(mPf 

c(k)c M 



M 



V / sup \V 2 g(xi,ty +{l-t) Xj )\dy 

Af 

E 

M 

^ , y^ oS<i 4^i^-%-(i-t)x,f d2/ 

r i 

<7iy 



3=1,3^ 
M 



3=1^1" 
M 



sup 

D, 0<t<l 



E 



)rf o<t<i 47r|5c; - ty - (1 - i)x-, | 3> 
1 1 



47r|a;i 


-ty 


-(1 


-t)xj\ 






1 




47T \xi 


-ty 


-(1 


~t)x 3 \ 2 






1 





(mP) 2 J 1 /( 2m p)<\x l -y\<V3 V 47r l^ ~ v\ ^\ x i ~ V? 4tt|^ z — 
c(k)c M 



dy 



(mP) 2 
2c(fc)cM 
(mP) 2 

2c(k)cM 
M 2 / 3 



r + 1 + - ) dr 

r , 



l/(2raP) 

l + log(V3)-log(l/(2mP)) 
1 + log(V3) - log 



= O 



/logM 



(101) 



where Cm := ||p|| 00 7 3 Ma^ ? ^, c(k) is a constant depending on the wave number 
k. Here the estimates \y — Xj\ < v / 3/(2mP) for y £ Dj, and |x; — y\ < 2\xi — s\ 
for y e Dj, j ^ I, s = tXj + (1 — t)y, i G [0, 1], were used. Using estimates (|95|) . 
([97). ([99]) and (fTUTjl . one gets relation ([?9"j>. 

Lemma T3. 31 is proved. □ 



The following theorem is a consequence of Lemma 13.21 and Lemma 13.31 
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Theorem 3.4. Suppose that the assumptions of Lemma \3.2\ and Lemma 
hold. Then 

\\u M - u eM \\ CM = O (jfj^ + |i _ 7 3 Ma 2 ~^ as M oo, (102) 

where um and u £ ^m are defined in (|61|) and (|80p . respectively. 

To get the rate of convergence (|102[) we have assumed that p(x) € C 2 (D). 
If p(x) € C(D) then the rate given in Theorem 13.41 is no longer valid. The rate 
of \\um — Ue.M\\c M when p(x) G C(D) is given in the following theorem. 

Theorem 3.5. Let Assumption A) hold and p € C(D) satisfies 
\p(x)-p(y)\<u p (\x-y\), Vx,yeD, (103) 
where u p is the modulus of continuity of the function p(x) . Then 

\\u M ~u eM \\ C M=0^^ + \l~ 1 3 Ma^\+u Jp (l/M 1 /^ , (104) 

where um and u e _M are defined in (|6ip and (|80[) . respectively. 
Proof. We have 

\\u M - U e ,M\\cM < ||«Af - Gm||c m + II^M - Me,Af||c M - (105) 

Let us estimate ||um — um\\c m - From (l67|) we have 

\\u M - um\\c m < ci||T d u - T d . M u M \\c M i (106) 

where ci is defined in (|52|) . Using the similar steps given in (|68|) we get the 
following estimate for \\TdU — T ( i ! m u m\\c m '■ 

\\T d u - T dM u M \\ C M < 3 l'Pll°°M°° +h+ h (1Q7) 

2(m/ J ) /: 



where Ji and /2 are defined in (|69|) and ([70]) . respectively. It is shown in (|72[) 
that Ji = 0(1/M 2 / 3 ). Since € C(D), the steps f|7S |) -([T7 |) are no longer valid. 
The estimate of -?2 can be derived as follows. Since p € C(D) and m € C 1 (D), 
it follows from (|70j) that 

M . 

h < IblUHfulU max V / |g(x ; ,y)||a; ; - y|rfy 

. V^lbHooll^ U||00 /" I i / -, ,, m ^ /71J-l/3\ 

< " " — y — max / —, :dy = 0(1 (mP)) = 0(1 M L/i ). 

(108) 
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This together with (fTUT|) and h = (9(1/M 2 / 3 ) yield 

\\u m -u m \\ c m =0(1/M 1 / 3 ). (109) 
Let us find an estimate for \\um — u e ^M \\c M ■ From (|84|) we have 

\\UM ~ U e ^ M \\c M < Cl\\T e U e ,M - 7d,MUe,M||c M ) (HO) 



where ci is defined in (|52|) . By definitions (|48|) and (|82|). and use the triangle 
inequality, we get 



\T e u e _ _ M - T rLM u 



M ~ J-d.M^e.MWC 1 



<Jl + J 2 



where 



J i := max 

KKM 



and 

J 2 '■= max 

KKM 



D, 



l{xi,y)p(y)dyu e {xi)dy 



A I 



I (9{xi,Xj)pa mP (xj) - g(xi,y)p(y))u e (xj)dy 



D, 



(111) 



(112) 



(113) 



p amp : = p(i)( 7 mP«( 2 - K )/ 3 ) 3 . It is proved in (|90j) that Ji = 0(1/Af 2 / 3 ). The 
estimate of J2 is derived as follows. By the triangle inequality we obtain 



J2 < J2S + J2,2, 
where 



J2 1 := max 

KKM 



and 



M 



^2 / (9(xi,Xj) - g(xi,y))p amP (xj)u e (xj)dy 



(114) 



(115) 



J2.2 := max ||u e || C M 

KKM 



It is proved in f|101[) that 
/logM 



A/ 

J] / l5(^>y)lba mP (a:j) -p(y)\dy. 



(116) 



J2.1 = O 



U 2/3 



(117) 
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To estimate we apply the triangle inequality and get 

M 

J2.2 < IKHc-w max V / \g(xi,y)\\p amP (xj) - p( Xj )\dy 

M . 

; 1^.7 J D4 



u P r u max 



M 



(118) 



< | 7 3 Ma^ p K - IIUpHoo max ||u e || C M V / \g(x h y)\dy 

3 = 1,3^1 J 
M . 

+ max sup u) p (\xj - y\)\\u e \\ C M max V / |p(x;, 
yen, ^^JD, 

= o(\ 1 3 Ma^-l\+u> p (l/M 1 / 3 )), 

where w p is the modulus of continuity of p(x). This together with 3\ = 
0(1/M 2 / 3 ) and (HHJ) yield 

\\um - u eM \\ C m = (j^ + | 7 3 Ma^ - 1| +c p (l/M 1 /3)^ . (119) 



Relation (pXi4]) follows from fl05| . (fT09"l) and (fTT9]) . 
Theorem 13.51 is proved. 

□ 

As we mentioned in the introduction, the main goal is to develop an al- 
gorithm for obtaining the minimal number of the embedded small balls which 
generate a material whose refraction coefficient differs from the desired one by 
not more than a desired small quantity. Let us derive an approximation of the 
desired refraction coefficient n 2 (x) generated by the embedded small balls. We 
rewrite the sum in (1381) as 



M 

^2g{x i x j )p amP {x j )u(x j )\D j \, (120) 
j=i 

where \x — Xj\ > a m p, j = 1,2,... ,M, \Dj\ = l/(mP) 3 is the volume of the 
cube Dj, and 



P, 



■a mP {x) := ATTh{x)N{x){ 1 mPa [2 - K ^ z f, N[x) = I/7 3 . (121) 



since as m — > 00, it follows that (|120[) is a Riemannian sum 

for the integral J D g{x, y)p{y)u{y)dy , where p(x) = 4irh(x)N(x). This motivates 
us to dehne the following approximation of the refraction coefficient n 2 (x): 

n 2 (x) := n 2 (x) - k~ 2 p amP {x), (122) 
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where p amP is defined in (I121[) . We are interested in finding the largest radius 
a m p (or the smallest M — (mF) 3 ) such that 

e(M) := max^ \n 2 ( Xl ) - n 2 amP ( Xl )\ < e/k 2 := e(k), (123) 

where k is the wave number, e > is a given small quantity and n 2 mP (x) is 
defined in (HH1). 

An estimate of the error e(M), defined in (|123[) . is given in the following 
theorem. 

Theorem 3.6. Suppose Assumption A) holds and N(x) — I/7 3 . Then 

max \n 2 ( Xl ) - I < fc- a |b|U 1 - hmPa^ /3 f , (124) 

l<K(mP) J 

where xi is the center of the l-th small ball, p{x) is defined in (|41[) , n 2 {x) = 
Uq(x) — k~ 2 p{x), and n 2 mP is defined in (| 1 22[) . Consequently, 

lim max \n 2 ( Xl ) - n 2 amP ( Xl )\ = 0. (125) 

m-s-oo l<l<( m P) 3 

Proof. Let 

h ■= \n(xi) -n 2 amp (xi)\. 
Then 

h = k- 2 \p( Xl ) - p(x0hmPa% p K)/3 } 3 \ < fc- 2 b(x ; )||l - [imPa {2 ' p K)/3 f\ 

(2-h 
1 mP 



<k- 2 \\p\Ul~l im Pc (2 - K)/3 ^ 



This together with relation (J22J) yield ([125)) . 
Theorem 13.61 is proved. 



(126) 



□ 

Using Theorem 13.61 one can calculate the smallest M satisfying (|123|) by the 
following algorithm: 

Algorithm 

Initializations: Let the wave number k, the constant e > 0, n§(x) and n 2 (x) 
be given. Fix P > 1, m = m := 1, k € (0,1), 7 > [1/(2P)]( K+1 )/ 3 and 
jV(x) = I/7 3 . Partition D into P 3 cubes ft g , D = U^J\, Slj n = for 
j i, where each cube flj has side length 1/P. 

Step 1. Solve the equation 

ia {2 p K)/3 + a mP - l/{mP) = (127) 
for a mP . 
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Step 2. Embed m 3 small balls of radius a m p in each cube Vl q so that Assumption 
A) holds. 

Step 3. Compute 

p(xi) = k 2 (nl(xi) - n 2 (xi)) 
and 

Va m M) = PMlimPai 2 n p K)/3 }\ 1 = 1,2,..., (mP) 3 , 

where x/ is the center of the Z-th small ball and k is the wave number. 

Step 4. If maxi<j<( m p)3 — Pa mP (£i)l > e ; then set m = 77i + 1 and go to 

Step 1. Otherwise the number M = (mP) 3 is the smallest number of 
the balls embedded in D such that inequality (|123[) holds, and a m p is the 
radius of each embedded ball. 

4 Numerical experiments 

In this section we give the results of the numerical experiments. Suppose the 
refraction coefficient of the original material in D is tIq(x) = 1 and the desired 
refraction coefficients are: 

Example 1. n 2 (x) = 5, 

Example 2. n 2 (x) — 5 + exp(— \x — x \ 2 / (2a 2 )) / where xq — (0.5, 0.5, 0.5) and 
a = Here b is the smallest number m taken from Example 1. 

Example 3. n 2 (x) = 1 + 0.5 sin(a;i), where x\ is the first element of the vector x, 

Example 4. n 2 (x) = 1 + 0.5 sin(100iEi), where X\ is the first element of the vector x. 

By the recipe we choose 

h{x) = k 2 [n 2 {x)-n 2 (x)]/(4irN(x)) = ~f 3 k 2 (n 2 {x) ~ n 2 (x)) / (Air) , k > 0. (128) 

Let us take 

P = ll, k = 0.99, 7 = 10/c[l/(2P)] (1+K)/3 , m = 1, (129) 

where k > 1 is the wave number and mo is the initial number of small balls 
described in the algorithm. Here the parameters P = 11 and mo = 1 are chosen 
so that the approximation error in Lemma 13.21 is at most c(fc)10~ 4 , where c(k) 
is a constant depending on the wave number k. We apply the algorithm given in 
Section 3 to get the minimal total number of small balls embedded in the cube D 
such that inequality f| 123[) holds for various values of e, where the quantity e was 
defined in the Algorithm (see the Initialization and Step 4 of the Algorithm). 
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The smallest number of the balls embedded in D increases as e decreases. 
The radius a rn p and the ratio a m p/d m p decrease as M increases, which agrees 
with the theory. The results are shown in tables 1-4. In these tables we define 

d mP := min dist(S QmP (x,), B amP (a;,)), (130) 
where B a (x) is dehned in (Ti~9")) . and 

E;= max \n 2 { Xl ) - n\ (a;,)!, (131) 

1<1<M a mP \ 

where M is the smallest total number of small balls embedded in the domain 
D, a m p is the radius of the embedded small balls and X[ is the center of the 
l-th small ball. In Example 1 we choose a constant refraction coefficient n 2 (x). 
For k = 1 the total number of small balls M increases by 1.651 x 10 5 when the 
error level e is decreased by 50%, while for k = 5 the value of M increases by 
3.4609 x 10 4 as the error level e decreases by 50%, as shown in Table 1. 

In Example 2 we add a Gaussian function to the constant refraction coef- 
ficient n 2 (x) considered in Example 1. For k = 1 the value of M increases 
significantly as the error level e decreases by 50%. 



k = 1 


e 






m 


M 


a„ 


iP 


a m p/d mP 
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5.000 x 


10" 


-1 


1 


1.331 x 10 3 


3.722 ; 


< 10" 4 


5.445 x 10~ 2 


9.7i7 x 


10" 


-2 


5.000 x 


10" 


3 


5 


1.664 x 10 5 
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< 10" 6 


1.098 x 10~ 2 


4.219 x 


10" 
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m 
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a m p/d m p 
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5.000 x 


10" 


-1 
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1.331 x 10 3 


3.200 : 


< 10" 6 


2.196 x 10~ 3 


8.448 x 


10" 
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5.000 x 


10" 


■3 


3 


3.594 x 10 4 
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< 10~ 7 


7.320 x 10 -4 


9.701 x 


10" 


-5 










Table 1: 


Example 1 










k = 1 
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1.331 x 10 3 


3.722 
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2.209 x 
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2.924 x 10 6 
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4.223 X 10~ 3 


4.715 x 
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k = 5 
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nP 


a mP /d mP 


E 






5.000 x 


10" 


-1 


1 


1.331 x 10 3 


3.200 


x ur B 


2.196 x 10~ 3 


1.915 x 


10" 


-3 


5.000 x 


10" 


■3 


5 


1.664 x 10 B 


2.686 


x 10" s 


4.392 x 10 -4 


1.702 x 


10" 


-4 



Table 2: Example 2 



The refraction coefficients n 2 (x) considered in Examples 3 and 4 are periodic. 
In Example 3 for k = 1 the total number of the embedded small particles M 
increases by 9319 as the error level e decreases by 50%. A similar increment of 
M is obtained for k = 5. These results are shown in Table 3. 

In Example 4 the angular frequency of the sine function is 100 times the 
angular frequency of the sine function given in Example 3. In this case we get 
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significant increments of the value of M as the error level e decreases by 50% 
for the wave numbers k = 1 and 5, see Table 4. 
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3.200 x 10~ 6 


2.196 x 10~ 3 


4.798 x 10" 


-6 


5.000 x 


10" 


-5 


2 


1.065 x 10 4 


4.084 x 10" 7 


1.098 x 10~ 3 


6.126 x 10 


-7 










Table 3: 


Example 3 








k = 1 


e 






m 


M 




&mP I dmP 


E 




5.000 x 


10" 


-1 


1 


1.331 x 10 3 


3.722 x 10~ 4 


5.445 x 10~ 2 


1.218 x 10" 


-2 


5.000 x 


10" 


-4 


6 


2.875 x 10 5 


1.861 x 10 -6 


9.147 x 10 -3 


3.673 x 10 


-4 


fe = 5 


e 






m 


M 


a m p 


a m p/d m p 


E 




5.000 X 


lo- 


-1 


1 


1.331 x 10 3 


3.200 x 10" 6 


2.196 x 10" 3 


1.05 X 10" 


i 


5.000 X 


ur 


-5 


8 


6.815 x 10 5 


6.650 x 10" 9 


2.745 x 10~ 4 


1.756 x 10 


-6 



Table 4: Example 4 
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